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Abstract 

The problem of the time required for a diffusing molecule, within a large bounded domain, to 
first locate a small target is prevalent in biological modeling. Here we study this problem for a 
small spherical target. We develop uniform in time asymptotic expansions in the target radius of 
the solution to the corresponding diffusion equation. Our approach is based on combining short- 
time expansions using pseudo-potential approximations with long-time expansions based on first 
eigenvalue and eigenfunction approximations. These expansions allow the calculation of correspond- 
ing expansions of the first passage time density for the diffusing molecule to find the target. We 
demonstrate the accuracy of our method in approximating the first passage time density and re- 
lated statistics for the spherically symmetric problem where the domain is a large concentric sphere 
about a small target centered at the origin. 
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I. INTRODUCTION 



Diffusion of a molecule to a spherical trap is a classical problem important in chemical 
kinetics. In an unbounded domain, the problem reduces to the Smoluchowski theory of re- 
action kinetics. In the context of biological processes, intracellular transport of biomolecules 
and chemical reactions occur within closed domains with complex geometries [6]. As a first 
passage time problem, this is closely related to the narrow escape problem, where a diffusing 
molecule escapes a closed domain through a small opening on the boundary, and the long 
time behavior has been studied using matched asymptotics [8, 10-12, 26, 27, 29]. There 
are many examples of this type of first passage time problem in biological modeling, in- 
cluding transport of receptors on the plasma membrane of a dendrite [4, 18], intracellular 
virus trafficking [23], molecular motor transport [5], binding of a transcription factor to a 
segment of DNA within a nucleus [21], and export of newly transcribed mRNA through 
nuclear pores [13]. 

Consider a bounded domain Q C 1R 3 , containing a small, absorbing spherical trap, Q e C Q, 
with radius e centered at rb £ fi. We denote by dQ the exterior boundary surface to 
Q, and by dQ e the exterior boundary to Q e . The non-trap portion of Q is denoted by 
f2f ree = O \ {Q e U dQ e }. Consider a molecule undergoing Brownian motion within f2f ree . We 
denote by p(r,t) the probability density that the molecule is at position r G fl{ ree at time t 
and has not yet encountered the trap. For D the diffusion constant of the molecule, p(r,t) 
satisfies the diffusion equation 

% = DV 2 p(r, t), re fi fl ,e, t > 0, (I.la) 
at 

d v p(r,t)=0, redtt,t>0, (Lib) 

p(r,t)=0, redQ e ,t>0, (Lie) 

p(r,0) = 5(r - r ), r G fi free , r G fi free (Lid) 

where d v denotes the partial derivative in the outward normal direction, 77, to the boundary. 
Let T label the random variable for the time at which the molecule first reaches dQ e . The 
first passage time distribution is defined as 



F(t) = Prob[T < t) = 1 - / p(r, t) dr. (1.2) 

Jn 

Suppose that the solution to (1.1) can be written in terms of an eigenfunction expansion 



with 

oo 

p(r, t) = J2 Mr )Mr)e- x "\ (1.3) 

n=0 

where the eigenfunctions and eigenvalues satisfy 

D\7 2 4> n (r) = X n ip n , r e O fr ee, (L4a) 

d v ip n = 0, re dn, (I.4b) 

Mr) = 0, re dn e . (I.4c) 

We order the eigenvalues so that < Ao < Ai < . . . . In the limit that the radius of the trap 
vanishes, the smallest eigenvalue, subsequently called the principal eigenvalue, also vanishes 
(i.e., Ao — > 0). Similarly, the corresponding eigenfunction, subsequently called the principal 
eigenfunction, approaches ipo(, r ) ~^ as e — )■ 0. Corresponding to these limits, the first 



passage time T — > oo and lim^oo f p(r,t)dr = 1 as e — > 0. In what follows we let 
diamS* and \S\ denote the diameter and volume of the set S C K 3 . For < e diamf2 
the asymptotics of the principal eigenvalue are known, and given by Ao ~ e Note 
that to first order in e, Ao depends only on the volume of Q and not the domain geometry. 
Higher order terms which depend on other properties of the domain are discussed in the 
next section. 

The small e asymptotics of Ao motivate a large-time approximation of p(r, t), based on a 
separation of time scales. Truncating the eigenfunction expansion (1.3) after the first term 
gives the long time approximation, 

p(r,t) ~ -pj e ~^ £ *> Ait > 1, e<diamfi. (1.5) 

Note, however, that the initial condition (Lid) is not satisfied by this expansion. Instead, 
the initial condition is modified so that the molecule starts from a uniformly distributed 
initial position with 

PM) = p- (1-6) 
In other words, the long time behavior depends very little on the initial position of the 
molecule because it is likely to explore a large portion of the domain before locating the 
trap. The long-time, Ait ^> 1, approximation of the first passage time density is then 

/(*) = j t T{t) ~ A e- Aot , Ait » 1, (I.7a) 

An De 4*p ct 
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e mi , e<diamO, (1. 7b) 
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where J-(t) is given by (1.2). The first passage time is therefore approximately an exponential 
random variable, with mean 



E[T] 




1 

V 
\n 



e <C diamfi. 



Xit > 1, 



(I.8b) 



(I.8a) 
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An exponentially-distributed first passage time is an important assumption in course-grained 
models, such as the reaction-diffusion master equation (RDME) [16, 17, 28]. (The RDME is 
a lattice stochastic reaction diffusion model which assumes that reacting chemicals are well 
mixed within a computational voxel.) More broadly, exponential waiting times are essential 
for jump processes to be Markovian. 

The above long-time approximation motivates several questions. First, when is the non- 
exponential, short-time behavior of the first passage time important? Second, how does 
changing the initial position of the molecule effect the approximation? It follows from (1. 8a) 
that the mean binding time is approximately independent of the initial position. On the 
other hand, as we show here, the most likely binding time, called the mode, depends strongly 
on the initial position. Recently, others have found that the initial position is also important 
when comparing two or more identically distributed binding times [24, 25]. More generally, 
to estimate spatial statistics for the position of the diffusing molecule it is necessary to 
obtain expansions of not just the first passage time density, f(t), but also the solution to 
the diffusion equation, p(r,t). 

Motivated by these and other examples, we develop a uniform in time asymptotic approx- 
imation as e of the probability density, p(r,t) (see (11.37)), and the first passage time 
density, f(t) (see (11.42)), that accounts for non-exponential, small time behavior and the 
initial position of the molecule. The paper is organized as follows. In Section II we further 
develop the long time approximation and present the complimentary short time contribu- 
tion based on a pseudopotential approximation. Combining these two estimates we derive 
an asymptotic expansion of p(r, t) for small e. In Section II C we use the results of Section II 
to derive a small e expansion of the first passage time density (through terms of order 0(e 2 )). 
Finally, in Section III these approximations are compared to the exact solution, exact first 
passage time density, and several other statistics for a spherical trap concentric to a spherical 
domain. 
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II. UNIFORM ASYMPTOTIC APPROXIMATION 



Our basic approach is to compute two asymptotic approximations for e <C diamfi, that 
are valid on two different timescales; "small times" Ait <C f, and large times Ait ^> f. We 
write p(r, t) as 

p(r,t) =p LT (r,t)+p ST (r,t), (III) 

where p hT is accurate for large times and p ST is the small time contribution that accounts 
for the initial position. If p LT is to be the large time approximation (1.5) with the initial 
condition (1.6), then the small time contribution must satisfy the projected initial condition 

p hT (r, 0) = (V>(r), 5(r - r )> V(r) = V^M^o), (IL2) 
p ST (r, 0) = 5(r - r„) - ^(r)V(r ). (IL3) 

Here we have dropped the subscript and subsequently identify ip and A as the principal 
eigenfunction and eigenvalue respectively. Using (II. 3) as initial conditions and setting t — 
in (II. 1) then gives us p(r,t) = 5(r — vq) as required. 

In the next two sections we derive asymptotic expansions of p hT and p ST for e <C diamfL 
The expansion of p LT is based off the principal eigenvalue and eigenfunction expansions 
developed in [9]. The expansion of p ST adapts the pseudo-potential method we first used 
in [20], where uniform in time expansions of p(r,t) and the first passage time distribution, 
Prob [T < t], were obtained for Q = IR 3 and the origin. We have found that a direct 
pseudo-potential approximation of (1.1) in bounded domains with Neumann boundary con- 
ditions breaks down for large, but finite times. For example, the direct pseudo-potential 
based expansion of Prob [T < t] can become negative for large times. For this reason we 
have introduced the long and short time splitting approach we develop in this work. 

A. Large time asymptotic expansion 

Since the initial condition (II. 2) is an eigenfunction of the Laplacian, the long time density 
is given by 

Pvr(r,t) = ip(r)ij(ro)e~ xt . 

As discussed in the Introduction, there are well known asymptotic approximations for small 
e of ip(r) and A. These then determine the small e behavior of p LT (r,t). The expansions of 
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i/j(r) and A are typically given in terms of the corresponding no-trap problem where e = 0. 
Let G(r, r', t) denote the fundamental solution to the diffusion equation in Q (i.e. the e = 
problem), then 





dt 



G{r, r', t) = DV 2 G{r, r', t), r 6 Q, 
G(r,r',0) = 5(r-r') 7 r e Q, 



(II.4) 



with the no-flux Neumann boundary condition 

d v G(r,r',t) = 0, re dVl. 



(II.5) 



We will also need the corresponding solution to the time-independent problem, the pseudo- 
Green's function U(r,r'), satisfying 

1 



DV 2 U(r,r') = -^--5(r -r'), rGd, 



(II.6) 



with the no-flux Neumann boundary condition 

d v U(r, r') = 0, re dVl. 
Note, by integrating (II. 4) we find that 

r ( G(r ' r/ ' t)_ pj) dt=u ^^)- 

Here the term is necessary to guarantee convergence of the integral. Finally, we denote 
by 7 the value of the regular part of U(r, r^,) at r = 77,, 



(II.7) 



7 = lim 

r— 



U(r,r h ) 



1 



(II.9) 



AnD \r — r^\ 

Let k = AtcD. As derived in [9], the asymptotic expansions of the principal eigenvalue 
and eigenfunction for small e are 

k 



\ ^ \ LT = — [I - k'je) e, 



(11.10) 



and 



V>(r) ^= + e^ (1) (r) + e 2 ij (2 \r), 



1 - ekU(r, r h ) - e 2 k 2 [ - 7 C/(r, r b ) + p J U(r, r')U{r', r h ) dr' 



(11.11) 



Here \1/ denotes the spatial average of the second order term and is given by [9] 

k 2 



(U(r,r h )Y dr. 



(11.12) 



2 Jn 

Note, the second order term in (II. 11) is not explicitly derived in [9], but can be found by 
solving equation (2.20) (of [9]) with the normalization (11.12). The corresponding expansion 
of the initial condition, Puvir, 0), in e is 

P LT (r,0) ~ p| + w {1 \r,r )e + w {2) (r,r )e 2 , 

where the functions w( n \r,ro) are obtained from substituting the asymptotic approxima- 
tions of the principal eigenfunction and eigenvalue into (II. 2) and collecting terms in e. We 
find that 



w (1 \r,r ) 



mi 



U(r,r h ) + U(r ,r b ] 



(11.13) 



w{2) ( r > r °) = iTTT ^( r > r b )U(r , r h ) + [U(r, r b ) + U(r , r h )} 



k 



k 2 ~t 



\n 



n\ 



\n\ 



f 2\I> 

/ [U(r,r') + U(r ,r')]U(r',r h )dr' + — =, (11.14) 

Jn V" 



so that the small e expansion of p^ir, t) is then 

1 



Pur(r,t) 



w^\r, r ) e + w (2) (r, r ) e 



(2), 



(11.15) 



B. Short time asymptotic expansion 

To construct an asymptotic approximation to p ST (r,t) for small e, we assume that the 
trap boundary condition can be replaced by a sink term in the PDE involving a Fermi 
pseudopotential operator [15, 19]. The boundary condition p ST (r,t) = for r G dQ € is 
replaced by the sink term 

d 



— ek 5(r — r b ) 



\r - r h \p ST (r,t) 



d \r — r h \ l 

For r = \r\, in the special case that r b = is the origin, this reduces to 

d 



(11.16) 



(11.17) 



The pseudopotential operator was developed for approximating hard core potentials in quan- 
tum mechanical scattering problems [15, 19]. In the case that Q = M. 3 , it has been shown 
that the asymptotic expansion for small e of the solution to the diffusion equation with 
pseudo-potential interaction approximates the solution to the diffusion equation with a zero 
Dirichlet boundary condition, p(r,t) = for r £ dQ e , up through terms of order 0(e 2 ) [20]. 
Here, we interpret the pseudopotential as an extension of the delta function to act on func- 
tions that have a singularity of the form r~ a for a < 1 (see [1-3, 14] for discussions of several 
approaches for rigorously defining pseudo-potential interactions). 
The pseudopotential approximation for p ST (r,t) is that 

d d r i 

— p ST {r,t) = DV 2 p ST (r,t) - ekS(r - r b )— , \r - r h \p ST (r,t) , (11.18) 

at a \r — rb| L J 

for r £ Q, with the initial condition (II. 3) and a no-flux Neumann boundary condition on 
dfl. Following [1, 14, 20] we split p ST (r,t) into a regular part, <p(r,t), and a singular part, 
q(t)U (r, r b ), so that 

p ST (r,t)=(f)(r,t) + q(t)U(r,r h ). (11.19) 

Here it is assumed that <f)(r,t) is "nice" as r — > r^. 

To find the asymptotic expansion of p ST (r,t) for small e, we begin by formulating a 
closed integral equation for <ft(r,t). Substituting the preceding representation of p 3T (r,t) 
into (11.18), we find 

^ = DVV - ^U(r, r b ) + q(t) (±- _ S (r - r b )) - efc (0(r b , t) + 7 q(t)) S(r - r b ). 

(11.20) 

We enforce the point boundary condition that the delta function terms should cancel [14, 20] 
so that 

g(t) = --^-0(r b ,t). (11.21) 
1 + e£ry 



(11.20) then simplifies to 



Using Duhamel's principle we find that 
4>(r,t) 



J Jj3{r,r',t-s) L^U(r',r h ) + j^q(sfj dr' ds 

+ [ G{r,r',t)(f){r',0)dr'. (11.23) 
Jn 



8 



Integrating by parts we find 

/ G(r,r',t-s) ( ^-ds= f DV 2 G(r,r',t-s)q(s)ds 
Jo ds J 



+ q(t)6(r - r') - q(0)G(r, r', t). (11.24) 



The no-flux boundary condition implies 

/ DV 2 G(r,r',t- s)U(r',r h )dr' 
Jn 

= [ G{r,r',t - s)DV 2 U(r',r h ) dr' 
Jn 

= jf G{r, r', t-s)Q^- 8(r' - r^j dr' 



mi 



G(r,r h ,t- s). 



(11.25) 
(11.26) 
(11.27) 
(11.28) 



Using the two preceding identities, it follows that (11.23) simplifies to 

4>(r,t) = -q(t)U(r,r h ) + / G(r, r h , t - s)q(s) ds 

Jo 

+ / G(r,r',t)[<P(r',0) + q(0)U(r',r h )]dr'. (11.29) 
Jn 

Eliminating q(t) with the point boundary condition (11.23) gives 
ke 



1 + 7/ce 



U(r,r h )<j>(r h ,t) - / G(r, r h , t - s)<j)(r h , s) ds 



+ / G(r,r',t)p ST (r',0) dr'. (11.30) 
Jn 



We now use the integral equation (11.30) to find an asymptotic expansion of p ST (r,t) in 
e. Let 

p ST (r, t) ~ pg?(r, t) + p«(r, t) e + pg?(r, t) e 2 . 
Similarly, we define the expansion of the regular part of p ST (r,t) by 

<f>(r,t) ~ {o) (r%t) + (1) M)e + (2) M) e 2 - 
Using (11.21) we identify the expansion terms of p ST (r,t) as 
p^(r,t) = ^°\r,t), 



p§2(r,t) = <P {1 \r,t) - k^°\r h ,t)U(r,r h ), 



„(2), 
PST 



(r,t) = <p {2 \r,t) + k f-0W(r b ,t) + kj^ \r h ,t)) U(r,r h ) 



The principal eigenvalue and eigenfunction expansions of the previous section imply that 

p ST (r, 0) ~ S(r ~ r o) - ~ w (1 \r, r ) e - w (2 \r, r ) e 2 . 
Substituting this expansion into (11.30) yields 

<f>(°)(r,t) = G(r,r ,t)--±- (11.31) 



4>^\r,t) = kU(r,r h )^°\r h ,t) - k / G(r,r h ,t-s)4> (0 \r h ,s)ds 

,/ ° (11.32) 

and 

0(2) (r } t) = -k 2 ^U(r, r b )0(°) (r b , t) + P 7 f G(r, r b , t - s)0 (o) (r b , s) ds 

Jo 

+kU(r,r h )4>W{r h ,t)-k [ G(r,r h ,t - s)^\r h , s) ds (11.33) 

- / G{r,r',t)w^{r',r )dr'. 
Jn 

Evaluating (11.33) requires the calculation of 0W(r b , t) = lim r _j. rb <p^(r, t). Let (jo(p, £) = 
G(r,r',t) — pr. Using (II. 8) we have that 

f/(r,r b )0(°)(r b ,t)- / G(r, r b ,t - s)^(r b , s) = 
J o 

[0(°)(r b , * - s) - 0^(r b , *)] G (r, r b , s) ds 

/oo i rt 

G (r,r h ,s)ds--^-J ^°\r h ,s)ds. 

Combining this expression with the identity 

r poo 

/ G(r,r',t)U(r',r h )dr' = / G (r,r h ,s)ds 
Jn Jt 

and reusing (II. 8) we find 

= - T [^°\r h ,t-s) -0 (o) (r b ,t)] G (r,r b)S )rf S 
A; Jo 

+ W\ Jt {O) ( rh > s " >ds + G ( rh ' r °' t " > J t G o( r ' r b, s)ds. 
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As r — > r b we then obtain 



r b ,t) = —k 



o 



G(r b , r , t- s) - G(r b , r , t) 
k 

+ m 



G (r h ,r h ,s) ds 

/GO /»00 
0< o Hr b ,s)ds + fcG(r b ,r o ,t) y G (r b , r b , s) ds. (11.34) 

Note, in the first integral Gq(t^, r b , s) will scale like s -3 / 2 as s — > 0. This singularity is 
weakened by the G(r b ,r ,£ — s) — G(r b ,r ,£) term, which formally scales like s as s — > 
(for fixed t > 0). As such, the overall singularity is s integrable. Similarly, G(r b ,r ,t) will 
cancel the effective singularity in t of the last integral. 
We therefore find the recursive expansion formula that 

Theorem ILL The asymptotic expansion of p ST (r,t) for e <C diamfi is given by 

pW{r,t) = G(r,r ,t)-j^, (II.35a) 

PST 



(r, t) = —k J G{r, r b , t - s)<j>V> (r b , s) ds + £v(r , r b ) (II.35b) 

I I «/ a£ 

pf T {r,t) = P 7 [ t G(r,r h ,t-s)^°\r h ,s)ds (II.35c) 
Jo 

-k ( G{r,r h ,t-s)(j> m {r h ,s)ds- [ G(r,r',t)w {2 \r',r )dr'. 
Jo Jn 

Let 

U = U(r ,r h ), f ST (t) = <P m (r h ,t) = G(r h ,r ,t) - p. (11.36) 

Combining Theorem II. 1 with the long time expansion (11.15) we find, 
Theorem II.2. For e <C diamfi, 

p(r,t) ~ G(r,r ,t) - (l - e^) - ^ [U(r,r h ) e - X ^ - (1 - e- A ^)f/] 

+ e W\ I G ^y^ t )W^)dr f -dz f G(r,r h ,t-s)f ST (s)ds. (11.37) 

Note, based on Theorem II. 1 one can derive an expansion of p(r,t) valid through terms 
of 0(e 2 ). That said, this expression is of sufficient complexity that we do not summarize it 
here. 
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C. First passage time distribution 



Denote by T the first passage time (FPT) for the diffusing molecule to exit through dQ t . 
The FPT distribution is defined as 



JF{t) = Prob[T < t] = 1 - / p(r,t)dr. 



(11.38) 



Substituting (II.l) into (11.38), we find that T{t) = 1 - ^/\tt\ip(r )e- xt - J n p ST (r,t)dr, 
where ip and A are the principal eigenfunction and eigenvalue satisfying (1.4) for n — 0. 
From (11.19) it follows that J Q p ST (r,t)dr = J n <f)(r,t)dr, so that 



T(t) = 1 - y/\n\if)(r )e 



-At 



(r, t)dr. 



(11.39) 



Denote by -7 r LT (t) = 1 — e Alt< the FPT distribution corresponding to the asymptotic ex- 
pansion of the long time approximation 1 — e~ xt (see Introduction and (11.10)). Notice from 
(11.31) that at leading order, f a <f>[r,t)dr ~ J n (p {0) (r,t)dr = 0. Substituting (11.11), (11.32), 
and (11.33) into (11.39) and collecting terms in powers of e yields 



1 - ekU - e z k 



2 1,2 



n\ 



U{r , r')U(r', r h )dr' - jU ) + 2e 2 ^ Vl^ 



ft i-t 

.2£,2. \ / t f„\ j„ , Jit I J,(l), 



+ [ek-e 2 k 2 ^J J f ST (s) ds + e 2 k j w (r b , s) ds, (11.40) 

where U and f ST (t) are defined in (11.36) and (j)^(r^,t) is given by (11.34). In evaluating 
the various spatial integrals we have made use of the identities J a G(r,r' ,t)dr = 1 and 
J n U(r,r')dr = 0. 

The FPT density function is then f(t) = 7j+J~{t). We denote the expansion of the long 
time scale approximation, Xe~ xt , by 



fw{t) = jfvr{t) = A LT e- A -* = A (l - he) ee-M'-H* 
(see (11.10)). Formally differentiating the asymptotic expansion of F{t) we find 

Theorem II. 3. The asymptotic expansion of f(t) for e <C diamf2 is given by 

1 



(11.41) 



/(*) 



I - eld 7 - rlr ( — J U(r , r')U(r', r h )dr' - jU ) + 2e 2 ^ VM 



hr{t) 



+ (ek-e 2 k 2 ^ f ST (t)+e 2 k(f>W(r h ,t). (11.42) 
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Since we have derived the expansion of f(t) by formal differentiation of the expansion of 
J-'(t), we obtain terms that are of higher order than 0(e 2 ) in (11.42) (as f LT (t) is 0(e)). In 
what follows, when referring to the "leading order", "first order", or "second order" expansion 
of f(t), we mean those terms arising from the derivative of the corresponding order expansion 
of J~(t), when treating J r LT (t) as 0(1). As such, the "leading order" expansion of f(t) will 
be fur(t), the "first order" expansion will be 

(l-ekU^)f LT (t) + ekf ST (t), 

and the "second order" expansion will be (11.42). 



III. A SPHERICAL TRAP CONCENTRIC TO A SPHERICAL DOMAIN 



To illustrate our asymptotic results we consider the problem of a diffusing molecule 
searching for a small spherical trap of radius e centered at the origin. We assume the trap 
is contained within a larger, concentric, spherical domain with unit radius. As this problem 
is exactly solvable, we will use the exact solution formulae summarized in this section to 
study the accuracy of our asymptotic expansions from the preceding sections as both e and 
the number of expansion terms are varied. 

Denote by p(r, t) the spherically symmetric probability density for a diffusing molecule 
to be a distance r from the origin at time t. We assume the trap is centered at the origin, 
so that r*b = 0, and let r = \r\, r = \r \. For p(r, 0) = 5(r — r )/r 2 , we have that 



p(r, t) 



p(r, t)dS, 



(III.l) 



9Bi(0) 



where dBi(0) denotes the boundary of the unit sphere. 

The advantage of this geometry is that an exact solution to the diffusion equation (1.1) 
is known [7]. We find 

oo 

p(r, t) = J2 a n <Pn{ro)<Pn{r)e~ Xnt , (III.2) 

where 



n=l 



<Pn{r) 



sm(x/Agl - r)) 



cos 



(VAn(l-r)) 



e < r < 1, 



a n = (0 n (r)) 2 r 2 dr, and the eigenvalue A n is given implicitly by 



tan~ 





nir 



0. 



(III.3) 
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The corresponding first passage time density is 

'M ~ -It 

where 



r 1 2 00 

/ p(r, t)r 2 cfr = — > b n X 

Je u n=0 



(III.4) 




e — cos (1 — el 




£sin Ul-e)J± 



\ 



1 + sin ((r - e) 



(l-e)(l + ^)-l 



An ^ \ 

D 



. (111.5) 



J 



In the remainder of this section, we list the quantities necessary to compute the asymp- 
totic expansions of p(r,t) and f(t) for small e. Recalling that = 0, U is then given 
by [29] 



U = U(r ,O) = -^-(- + r 4- 9 -). (ITI.G) 



It follows from (II. 9) that 
and from (11.12) that 
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The fundamental solution G(r , 0, t) = g(ro, 0, t)/4n, where g(r, r , t) denotes the spherically- 
symmetric Green's function for the e = Neumann problem (see Appendix A), is given 
by 

1 00 -.00 
G(r , 0,t) = — + Yl ^ e ~^ °> *) = TnT + E a " e ^' 



n=l 



n=l 



where 



a„ = — [1 + — J = </„ sine 




r 



with sinc(x) = sin(x)/x. The eigenvalues, /x n , satisfy 

r"n \ / /^n 



tan 
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+ mr = 0. 
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Note that by comparing (III. 3) to (III. 10) it follows that lim^o = £*n- Integrating (11.37) 
over the unit sphere we find 



p(r, t) ~ g(r, r , t) - 3 (l - e" ALT *) - 3e£; [C/(r, 0)e" 
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The asymptotic expansion of the first passage time density, /(£), can be evaluated directly 
from (11.42). Here we use (II. 8) and (III. 8) to express U(r, r') as an eigenfunction expansion. 
We find that 



/ U(r ,r')U(r',0)dr' = T^. 
Jn n=1 t* n 

The short time contribution to the first passage time density is given by 

„ t oo 
Jo ti^n 

/oo / 1 °° \ oo 

G (r h , r b , s)ds = I -—- + ^e'^ 1 Yl 
\ ' ' n=l / m=l 



while 



To evaluate the time convolution 



G(r h , r , t- s) - G(r h , r , t) 



G (r h ,r h ,s)ds 



(111.12) 



(111.13) 



(111.14) 



in (11.34) we use the Python quad routine. The integral is split into a short time portion, 
s G (0, s*), and a long time portion, s G (s*, t). s* is chosen sufficently small that G(r^, s) 
can be approximated by a Gaussian evaluated at the origin, (4vr J Ds)~ 3/2 , with the same 
absolute error tolerance we use in evaluating the preceding series (see Appendix B). 



A. Results 



We now study the error between the exact spatial and first passage time densities from 
the preceding section, p(r,t) and f(t), and their asymptotic approximations for small e. In 
what follows we keep R — 1, D — 1, and vary e between 10 -4 and 10 _1 . The tolerances we 
used in evaluating the various series of the previous section are given in Appendix B. 

While we are interpreting our spatial and time units as non-dimensionalized, these choices 
are also consistent with using spatial units of /mi and time units of seconds. With these 
units the overall domain has roughly the radius of a yeast cell nucleus. We may therefore 
interpret the trap as a DNA binding site that a diffusing protein is searching for. While 
trap radii for DNA binding sites are not generally experimentally measured, the width of 
some DNA binding potentials have been measured. For example, the LexA protein binding 
potential was found to have a width of approximately .5nm [22]. 
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FIG. 1. Relative error in approximating the principal eigenvalue, A, by A LT . Observe that the error 
decreases like e 2 as expected from (11.10). 

The long time approximation of the first passage time density is the single exponential 
Aexp(— At), with the timescale A -1 . The principal eigenvalue A is given implicitly by (III. 3) 
(with n — 0) and has the asymptotic approximation A ~ A LT (see also (11.10)). Hence, 
for small e, the long time approximation of the first passage time density is asymptotic to 
fur{t) = A LT exp(— A LT t). As described at the end of Section II C, we refer to /lt(^) as the 
leading order approximation of f(t) as e — > 0. (We will also interchangeably refer to / L t(£) 
as either the large time or long time approximation.) 

The implicit equation (III. 3) can be solved numerically to calculate A to arbitrary preci- 
sion by a root finding algorithm (e.g., Newton's method). In Fig. 1 we compute the relative 
error, |A — A LT | A -1 , of the asymptotic approximation, A LT , as compared to the numerically 
estimated value of A (computed to machine precision). We see that as e — > 0, the relative 
error between the two decreases like e 2 , as expected from (11.10). 

In Fig. 2, we show the leading order spatial density approximation (blue curve), the first 
order expansion (green dashed curve), and the exact spatial density (black curve). These 
curves plot 

p(°)(r, t) = g{r, r , t) - 3 + 3e" ALT *, (111.15) 

the expansion (III. 11), and p(r,t) (III. 2) respectively. The spatial density is shown as a 
function of r at four different time points. For this figure we set r = 0.8 and e = 0.001. 
The density is initially concentrated around the initial position at t = 0.001 and slowly 
fills the region e < r < 1 until the density is approximately uniform at t — 1. The only 
visible difference between the leading order approximation and the exact result is near the 
absorbing boundary, r = e, where the exact solution displays a boundary layer that is lost 
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FIG. 2. The spatial density function p(r,t) (black curve) and its asymptotic expansions for small 
e at several time points. The blue curve gives the leading order expansion (III. 15), p(°)(r,t), while 
the green dashed curve gives the first order expansion (III. 11). We use a logarithmic x-axis to 
emphasize the solution behavior near the target, ro = 0.8 and e = 0.001 (similar to the width of 
measured DNA binding potentials [22]). 

in the leading order approximation. The first order expansion (III. 11) reintroduces this 
boundary layer and is indistinguishable from the exact solution at the scale of the graph. 

In the remainder of this section we focus on the approximation of the first passage time. 
The only free parameters in the model are the radius of the trap, e, and the initial distance 
from the trap, r . The long time approximation, f LT (t), is independent of r . It follows that 
the accuracy of fvr(t) in approximating f(t) improves when the initial distance from the 
trap is large (i.e., e <C ro < 1). In other words, the long time approximation is best when 
the particle is likely to explore a large portion of the domain before locating the trap. When 
the initial distance from the trap is small (i.e., e < r C 1), we might expect the short time 
contribution to be significant since there is a higher probability that the particle will quickly 
locate the trap before exploring the rest of the domain. In Fig. 3, we show the asymptotic 
expansion of the first passage time density (11.42) for e = 0.05 and ro = 0.3. With this choice 
the initial distance of the particle from the trap is small. Moreover, since the accuracy of 
the expansion (11.42) should decrease as e increases, taking e = 0.05 demonstrates the worst 
case behavior of the expansion for biologically relevant values of e. 

In Fig. 3 (left) the density function is shown with t on a log scale to accentuate the 
small time behavior. There is a significant difference between the long time approximation 
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FIG. 3. The first passage time density, f(t), for ro = 0.3 and e = 0.05. Asymptotic approximations 
of varying order are compared to the exact solution. The left plot uses a logarithmic i-axis and 
linear /-axis, while the right is linear in t and logarithmic in /. 
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FIG. 4. The first passage time density, f(t), for ro = 0.8 and e = 0.05. Asymptotic approximations 
of varying order are compared to the exact solution. See Fig. 3 (left panel) for the legend. 

(light blue curve) and the exact solution (green curve). The first and second-order uniform 
approximations correct for this difference. The large-time behavior is shown in Fig. 3 (right) 
with / on a log scale. For all except the shortest times the curve is linear, reflecting the 
exponential long-time behavior. We see that on this timescale there is very little visible 
difference between each curve. Fig. 4 is the same as Fig. 3, except that r = 0.8 so that the 
initial distance from the trap is larger. In this case the peak in the density occurs at a larger 
time. In both cases, the qualitative difference between the exact solution and the long- 
time approximation is a time lag before the exponential long time behavior dominates. The 
timescale for this time lag is roughly the diffusive transit time to cover the initial distance 
from the trap (i.e., r^/D). 

The absolute error of these approximations is shown in Fig. 5 for r = 0.3 and ro = 0.8. 
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FIG. 5. Absolute error of the first passage time density approximation for ro = 0.3 and vq = 0.8 
with e = 0.05. 
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FIG. 6. The max norm error of the approximation as a function of e. Solid curves show the 
error of the long-time approximation fur{t). The dashed curves show the second order uniform 
approximation. Note that the r$ = 0.65 and ro = .9 curves for the large-time approximation are 
indistinguishable . 

In both cases, the maximum error is noticeably decreased as the order of the asymptotic 
expansion is increased. Comparing the first and second order expansions, we see the main 
increase in accuracy results for times less than t — 1. Points in time where one of the 
approximations crosses the exact solution result in locally increased accuracy (the cusp-like 
drops in the expansion errors). Interestingly, when ro = 0.8 the long time approximation is 
more accurate for large times than the first- or second-order uniform approximations. Note, 
however, the error in each expansion at these times is substantially smaller than for short 
to moderate times. 

Finally, we examine the max norm error, max t > | /exact if) — f(t)\, as a function of e for 
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FIG. 7. The mode of the binding time distribution, defined as the most likely binding time, as a 
function of tq. Solid curves show the exact solution, dash dotted curves the first order approxima- 
tion, and the dashed curves the 2nd order approximation. 

different values of tq. The time points this error was numerically evaluated over are the same 
as those used for the graphs in Fig. 5, and are given in Appendix B. The result shown in 
Fig. 6 confirms the asymptotic convergence of the approximation as e — > 0. The large-time 
approximation (11.41) error (solid lines) shows linear convergence, while the second order 
uniform approximation (11.42) error (dashed line), which includes short time behavior, shows 
cubic convergence. 

As stated in the Introduction, the mean binding time is well approximated by the r - 
independent large-time approximation. That is, E[T] ~ 1/A, where A is given by (11.10). 
However, other statistics may be of interest that depend strongly on ro. One example is 
the mode, denned as the most likely binding time, call it r m , where f(r m ) = maxo<t<oo /(£)• 
Since the large-time approximation is an exponential distribution, the corresponding ap- 
proximation of the mode is r m ~ 0. In figure Fig. 7, we compute the mode by numerically 
maximizing the first passage time density. The exact mode is compared to first (dash dot- 
ted curves) and second order (dashed curves) approximations of the mode as a function of 
To- Each of the indicated curves are drawn for three different values of e. For e = 10~ 3 , 
the difference between each curve is indistinguishable. Notice that as e decreases the mode 
increases, particularly for larger values of Tq, indicating that the large time approximation 
of the mode becomes less accurate as e — > 0. 
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IV. DISCUSSION 



Although the first passage time of a Brownian particle in a confined geometry is a well- 
studied problem, an analytical characterization that includes short-time behavior of the 
survival probability density has been unresolved. The asymptotic approximation of the 
long-time behavior establishes a link between the spatial characteristics of the problem (i.e., 
the starting position of the particle and the space dependent survival probability density) and 
the short time behavior. That is, the long time approximation loses information about the 
initial position and treats the survival probability density as uniform in space. Consequently, 
the long time approximation is insufficient if one is interested in statistics that depend on 
these spatial characteristics. 

Using a multiple time-scale perturbation approach, we develop separate long and short 
time expansions of the solution to the diffusion equation in a bounded domain containing a 
small, absorbing spherical trap. The long time approximation is derived from the matched 
asymptotic expansions of [29], while the short time expansion is derived by modification 
of the pseudopotential method used in [20]. Combining these expansions we develop a 
uniformly accurate (in time) approximation of the survival probability density and the first 
passage time density. To study the accuracy of our method, we consider a example problem 
where the domain and trap are concentric spheres. By assuming radial symmetry, we have 
available for comparison the exact solution to the example problem. Our results show 
excellent quantitative agreement for all times over a range of physiologically realistic values 
of e. Moreover, they demonstrate the applicability of our expansions to estimating statistics 
that depend critically on the initial position of the diffusing particle. 
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Appendix A: Spherically-symmetric Neumann Green's function 



Let g(r,r ,t) denote the spherically symmetric solution to the diffusion equation, satis- 
fying 

2 d 9~ 



dt r 2 dr 
dg 



dr 



dr 



0. 



r e [0, 1) , 
r = 1, 



with the initial condition that g(r, ro, 0) = 5(r — ro)/r 2 . With this choice, 



g{r, r ,t) 



G(r,r ,t) dS, 



'aBi(o) 

for dBi(0) the boundary of the unit sphere. Here G(r,ro,t) denotes the solution to the 
corresponding three dimensional diffusion equation (II. 4). Note also the normalization that 

-i 



g(r, ro, t)r dr = 1 




G(r, r , t) dr. 



By eigenfunction expansion we find 



oo 

s(r,r ,t) = 3 + 2^(l 



n=l 



—J smc I W —r 1 smc 




where the eigenvalues /i n satisfy (III. 10) and we use the convention that 

• / x s in (a ; ) 
smc (x) = . 



x 



Appendix B: Numerics 

When evaluating the series for the exact solution (III.4) and asymptotic approxima- 
tion (11.42), we sum until the magnitude of the last added term drops below a given error 
threshold. We used an error threshold of 10 -14 for the exact solution and 10 -7 for the uniform 
approximation. The figures are generated with 1000 equally-spaced points for 10~ 3 < t < 1 
and 500 equally-spaced points for 1 < t < 30. 
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